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Abstract 

I propose a method, based on a set of Langevin equations, for bringing classical 
gauge theories to thermal equilibrium while respecting the set of Gauss' constraints 
exactly. The algorithm is described in detail for the SU(2) gauge theory with or 
without the Higgs doublet. As an example of application, canonical average of the 
maximal Lyapunov exponent is computed for the SU(2) Yang-Mills theory. 
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1 Introduction 



In studying thermal properties of gauge theories it is useful to consider the classical limit. 
This is particularly true for classical lattice gauge theories whose numerical simulation 
is increasingly used as a nonperturbative tool for real-time dynamics of the gauge fields 
(sphaleron transitions, properties of quark-gluon plasma 0, [3], f|, H || §). 

Simulation of a finite-temperature field theory should involve generating the ther- 
mal ensemble of field configurations. In case of gauge theories, configurations com- 
prising the ensemble are subject to a set of local Gauss' constraints, whose presence 
precludes straightforward application of standard importance-sampling methods of en- 
semble generation in most cases. This work deals with construction and implementation 
of constraint-respecting thermalization algorithms for a number of physically important 
classical gauge theories. 

So far, three methods have been used to thermalize gauge theories. One method 
consists of explicitly solving the constraints and applying importance sampling to the 
remaining gauge-invariant variables |J. Its practical value is probably limited to simple 
Abelian models. For more complicated, non-Abelian theories the dynamics in a gauge- 



invariant language is usually nonlocal and/or suffers from coordinate singularities |L6 
The second method handles Gauss' constraints approximately by adding to a Hamil- 
tonian terms which penalize deviations of the static charge density (SCD) from zero 
|2|, [|. Configurations with large SCD would then be unlikely to appear in a sample 
generated by an importance sampling procedure. A special cooling procedure is used to 
remove whatever SCD has nevertheless been generated in the system. This method has 
its shortcomings too. Large S CD-suppressing terms will dominate the energy functional 
and slow down the importance sampling. Besides, cooling may lead to deviations from 
the intended thermal ensemble. Yet another method circumvents the need for preparing 
initial thermal configurations by employing a self-consistent heat bath interacting with 
a finite-size system at the boundaries in real time [^, This method, like the first one, 
faces technical difficulties when applied to nonabelian theories or in dimensions higher 
than one. 

One would like to have a thermalization algorithm formulated in terms of usual 
phase-space variables of a gauge theory (e.g., gauge potentials and color electric fields) 
and accurately maintaining the Gauss' law. It turns out that such a method can in several 
important cases be based on Langevin equation with a specially designed multiplicative 
noise term. That Langevin equations with multiplicative noise can be used for thermal- 
ization of constrained systems has been known for quite a long time now. This is how 
unitarity of link variables is maintained in Langevin simulations of Euclidean gauge the- 
ories fK| . The unit-magnitude constraint on fields in a sigma model can also be treated 



this way [T^, O]. In what follows the same idea is employed in order to satisfy the set 
of Gauss' constraints in a gauge theory. As we shall see, a natural choice of algorithm 
in this case reflects the fact that Gauss' constraints are first class in Dirac's terminology 



|L1| . For technical reasons, it is convenient to use for our purpose second-order Langevin 



equations, i.e., of every pair of canonically conjugate phase-space variables only one will 
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be coupled to a heat bath by adding stochastic terms to its equation of motion [12 



The outlined strategy is discussed in some detail in the sections to follow. In Section 
2 I develop a Langevin formalism for Hamiltonian theories with first-class constraints. 
In Section 3 constraint-respecting coupling to a heat bath is constructed for a number of 
gauge theories. Based on this construction, a Langevin dynamics is formulated for the 
SU(2) lattice theory with and without a scalar doublet in Section 4. Section 5 discusses 
numerical integration of the resulting system of Langevin equations. Application of the 
Langevin algorithm is illustrated in Section 6 where canonical average of the maximal 
Lyapunov exponent is computed for the SU(2) Yang-Mills theory. Section 7 contains 
conclusions and outlook. A brief outline of the ideas presented here was given in Ref . Pf . 



2 Langevin equations with first-class constraints 

In order to present the method in the most general form, I must briefly review the basic 
phase-space properties of gauge theories. In doing so, I will closely follow Ref. p5| . 
Let a dynamical system be described by n coordinates g« and their respective conjugate 
momenta pi forming a phase space T. Let these variables be subject to m (m < n) 
time-independent first-class constraints 

C a (p,q) = 0^{H,C a }^{C^C a }, 

where H (p, q) is the Hamiltonian. The constraint surface M. is a 2n — m-dimensional 
subspace of the phase space. The weak equality sign, means, as usual, "equal up to 
terms vanishing on the constraint surface". A quantity f(p,q) is physically meaningful 
(observable) if it is gauge- invariant: {/, C a } ~ 0. Obviously, the Hamiltonian H is 
observable. Since any such / is constant along the gauge orbits generated by C a in A4, 
the physical subspace T* of the phase space is of dimension 2(n — m). Since the Poisson 
bracket of any two observables is again an observable, T* is by itself a phase space with 
a Poisson bracket {}* defined for observables u and v as 

{u,v}* = {u,v}\ M . 

It follows that there exists a canonical basis p*,q*,l < i < m for functions in T* such 
that 

M,P*}* = {<?*,<?*}* = 0; K, = 

All the observables (and H in particular) on M. depend exclusively on x* = p*,q*. 

Consider now Langevin dynamics in T. In its most general form the system of 
Langevin equations is 

x i = {H,x i } + D i {x)+gi(x)ri(t)), (1) 

where I collectively denoted p and q by x. The system (||) is intentionally written 
as a generalization of Hamiltonian equations of motion: the first term on the r.h.s. 
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describes the canonical evolution, whereas the other two terms arise from interaction 
of the dynamical system with a heat bath. The Di terms are, in general, dissipative. 
The last term on the r.h.s. is proportional to a white-noise random variable T J (t) with 
average zero and correlation 

(Ti(t)T k (t>)} = 5 jk 5(t-t>). (2) 

Here and in the following terms of this kind are to be understood in the Stratonovich 
sense, unless otherwise stated. At the moment, I do not specify the number of inde- 
pendent random variables T^(t). Our goal of thermalizing the gauge theory by means 
of will be achieved if (a) it preserves the constraints, i.e., generates a sequence of 
configurations on M. from an initial condition on Ai, and (b) for a long evolution time 
generates a sequence of configurations in T* distributed with canonical probability den- 
sity exp(—/3H(p*, q*), where (3 is the inverse temperature. Since T J are random variables 
independent of x, condition (a) requires that 

ddiC a « 0. (3) 

Now I will show that condition (b) can be satisfied together with (a) for a suitable choice 
of Di and g\. In order to arrrive at such an Ansatz, consider first the Fokker-Planck 



equation [18 for the probability density W(x,t) in the full phase space T, which follows 
from ([J) 

~d t + d t ({H, x t } + Di + g{d k gf - d k g{gi)} W(x, t) = (4) 

and require that it have exp(— f3H(x)) as its static solution. To this end it is enough to 
choose Di = g k G k , where 

G k = d j9 * - PgfaH. (5) 

With this form of Di (|l|) is constraint-preserving if (|3|) holds. Our task is now reduced 
to finding a correct set of en which in turn would determine Di. To further narrow 
down the search for a suitable I require that ([]]), averaged over realizations of random 
variables, maps observables to observables. The equation of motion for the average 
(f(x,t)) = J d 2n xW(x,t)f(x) of an observable f(x) is obtained from (|J) integrated by 
parts over the full phase space with boundary terms discarded: 

(/) = ({H, /} + d k (glgidj) - JglgphM'hf). (6) 

Gauge invariance of (/) is ensured by choosing 

gl = M^ k (x){T k (x),x t }, (7) 

where T k is observable and the square matrix M.M^ k is such that MM T is also an 
observable quantity. Indeed, with this choice the quantity averaged over on the r.h.s. of 

(§) i s 

{H, /} + {T k , M jk M jl {T l , /}} - f3M jk M jl {T k , H}{T l , /}. (8) 

In addition, gj as given by (0) satisfies @. If the initial configuration for (|l|) lies on 
/A, (PD, being an observable, depends exclusively on gauge- invariant canonical variables 



4 



x*. Moreover, Poisson brackets {} on T can be replaced by their counterparts {}* on 
T* everywhere in (§). Finally, gauge invariance of @ allows introduction of a gauge- 
invariant probability density W*(x*,t) on T* such that (/) = / d 2ijl ~ ni) x*W* f for any 
observable / on M.. In view of @ and (jD W* satisfies the Fokker-Planck equation 

8 t W* = {W*, H}* - {T k , M jk M jl ({T l , W*}* + (3W*{T\ #}*)}*, (9) 

written entirely in terms of x*. It is obvious that (^) has exp(— (3H), with H restricted 
to A4, as its static solution. In the following I shall reserve the terms "generators" for 
T k and "multiplier" for M jk . 

Substitution of (|7|) into ([!]) gives the equation of motion for an arbitrary variable v: 
v = {H, v} + ({T fc , M jk } - (3M jk {T k , H} + T j ) M jl {T l , v}. (10) 



Note that, even though ([I]) was written in a canonical basis, (|10D is basis- independent. 
This property, familiar from the Hamiltonian equations of motion, is useful whenever 
the natural choice of independent dynamical variables is not canonical, as is the case 
for Hamiltonian lattice gauge theories. Another observation is that it is always possible 
to include in (|10D terms with T k = C a . One one hand, this is allowed because the 
constraints are themselves first-class quantities. On the other hand, such an inclusion 
obviously does not change the equations of motion for observables and is therefore a 
matter of convenience. 

There are important special cases of (0). If M- 7 = 5^ k , the r.h.s. of (|T0| ) is observable 
for any observable v and separately for any realization of T^, not only on the average. In 
this case the stochastic terms in (|I],[lT]) m ay be thought of as arising from interaction of 
the dynamical system with gauge-invariant degrees of freedom of the heat bath. Some 
of the examples to follow are of this kind. 

The number N of independent random variables T- 7 and the functional form ofT fe , M^ k 
must be specified so as to optimize the convergence of W*(x*,t) to exp(— (3H). If N — 
2(n — m) it is always possible to choose M jfc = 5^ k , T k = x* k for any canonical basis x* 
in r*. With this choice the heat-bath fluctuations described by generate the most 
general variation of variables in T*. The system (|I|) then becomes a standard system 
of Langevin equations with additive noise on T*, known to guarantee convergence to 
the equilibrium distribution. Thus the class of Langevin equations defined by (|5|) and 
(0) is rich enough to achieve thermal equilibrium in any gauge theory. Convergent and 
numerically competitive algorithms are obtained by taking N = n — m, T k = pi, i.e., by 
coupling the heat bath to a single variable in every canonically conjugate pair p* , q* [|T2] . 



In the following these algorithms are called second order. In the next section I construct 
Langevin systems of this type for a number of gauge field theories in the continuum. 
Later on I discuss analogous algorithms for lattice gauge theories. 
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3 Langevin systems in the continuum 



In this Section I study how the general prescription for finding a correct heat bath 
coupling, as outlined in Section 2, works in a number of gauge theories in the continuum. 
It is convenient to set, with no loss of generality, A = 0. I will assume space dimension 3, 
which is the case for the most important applications. It is implied that all phase-space 
variables carry space coordinate labels, omitted here in order to simplify the notation. 
It should be clear that Langevin systems of equations for a classical field theory in 
the continuum are formal and cannot be put to immediate use: the Ray leigh- Jeans 
divergency renders the classical thermal theory in the continuum ill-defined. However, 
the continuum Langevin systems are useful as a guidance for similar constructions on a 
lattice. 

The simplest example is that of a free electromagnetic field. The Gauss' law involves 
exclusively the electric fields Ei and reads diEi = 0. Hence an arbitrary variation of 
the gauge fields A4 satisfies the Gauss' law. Such an arbitrary variation is generated by 
Ek, k = 1,2,3 and can be decomposed into a variation of the physical, transversal part 
of the gauge potential and a gauge transformation. If T k = Ek, M^ k = 5^ k in (|i0|) , it is 
straightforward to verify, by writing (|TI]) for A and E in momentum representation, that 
the resulting system of equations is indeed second-order for the transversely polarized 
degrees of freedom. 

Next, I discuss the SU(2) Yang-Mills theory. The Gauss' constraints for color gauge 
fields Af and momenta Ef are 

diE? - 2e a ^A^E] = 0, (11) 

where a, 0, 7 are adjoint color labels. For fixed electric fields the most general variation 
of the gauge fields can be written as 

6A? = SijEf, (12) 

where Sij is a (space-dependent) gauge-invariant tensor. Indeed, if the 3x3 matrix 
is invertible, as is the case generically, we find SV,- = dAfi^E' 1 )". For the Gauss' law ([TlD 
to hold, however, must be symmetric. Variations of A of this form are generated by 
the six quantities = EfE?, where [ij] labels ordered pairs with 1 < i < j < 3. The 
variables are gauge invariant, generically functionally independent, and commuting 
with each other. Hence they can be chosen as momentum variables in the observable 
subspace of the phase space. In fact, it is not difficult to find the corresponding gauge- 
invariant conjugate variables. To this end the original color electric fields are expressed 
in terms of T™ and three Euler angles 6 which parametrize an orthogonal matrix R aa {9) 
transforming E? to some standard form For instance, as pointed out by Goldstone 



and Jackiw [TH] , it is possible to require that be a product of an orthogonal matrix 



by a diagonal one. Then for any X = T^, 9, 

mi ^ d 



v , = A]^- - e^^- (R™(6)d t R^(6)Ej) 
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is the canonically conjugate variable. Corresponding to the described canonical trans- 
formation is the generating function 



$(A,X) = AjEj{X) - e afh R? a (8)diR Pa (6)E% 

whose partial derivatives give canonical conjugates of its arguments. Gauge invariance 
of Pt can be readily verified by performing a gauge transformation 

with an orthogonal Q. Consequently, the choice of = EfE^ as generators and 
= as a multiplier results in a second-order Langevin system in the physical 

subspace of the SU(2) Yang-Mills theory. The system of equations for the standard 
variables can A, E can be easily derived using (10) . It is obvious that the E equations 
are purely Hamiltonian. 

Inclusion of scalar fields in the scheme presents no difficulty The simplest theory of 
this kind is scalar electrodynamics. In this case, the Gauss' law is a relation between the 
electric field and the charge density of the complex scalar field 

C = diEi + i{ir<f>-<f>*Tr*) = Q, (13) 

where <fi and 7r are scalar field and momentum, respectively. The constraint C generates 
gauge transformations 

Ai — > Ai — didc] 4> — > 0exp(— ia); n ^ n exp(ia), 

where a is a gauge function. The most general variation of fields consistent with ( |T3"D 
consists of an arbitrary variation of Ai and a change in proportional to 7r*: 

5Ai = Vi\ 5<j) = S7i*. (14) 

These are generated by E i: i = 1, 2, 3 and II = tt*ti,. It is easy to identify the four (per 
space point) physical degrees of freedom varied by ([l4]) . Obviously, one can choose Ei 
and the magnitude II of 7r as gauge-invariant momentum variables. The fifth momentum 
variable, the phase 9 of 7r, is clearly gauge-dependent. Gauge-invariant canonical partners 
of Ei and IT follow from the generating function 

${A, 0, E, n, 9) = EiAi + vr(n, 9)<P + 0*7r*(n, 9) - EM- 

These are 

Ai — di9 and — -= (<fiexp(i9) + <fi* exp(—i9)) , 
2v n 

respectively. Choosing Ei, H as generators and a unit multiplier completes the construc- 
tion of a second-order system for the scalar electrodynamics. 

My final example, the SU(2) theory with a scalar doublet, describes a sector of the 
Standard Model. The Gauss' law now reads 

C a = &E? - le^A^Ej + i{na a <f) - 0V q tt*) = 0. (15) 
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Here the SU(2) spinors and 7r are scalar fields and momenta, a are Pauli matrices, and 
the rest of notation is obvious from previous examples. 

As before, in order to find a suitable set of generators and multipliers for this system, 
it is useful first to identify the physical degrees of freedom. To this end, a point transfor- 
mation of momenta is performed. Introducing II = tttt*, the momentum transformation, 
possible in the generic case II 7^ 0, reads 



71 = (VTl 0)lP{B); Ef = ^Tr {a a U{9)£^a^{9)) = R af3 {9)£?, 



where the Euler angles 9 parametrize the unitary matrix U rotating tt to a standard form, 
here chosen to be (IT 1 / 2 0). Matrix R, defined by the last equality, is the orthogonal 
representation of U. Only 9 change under gauge transformations, while II and the nine 
£ variables are gauge- invariant. As in the previous examples, the canonical conjugates 
of IT and £ obtained with the help of a generating function 



0, 0*; £, n, 9) = A?E?(£, 9)+7r(U, £)0+0*7r*(II, 9) + -Tr (£^^(9)8^(9)) (16) 

are again gauge invariant. These are 

n = (vr0 + 0*7r*)/2n 

for II and 

At = ^Tra a (u\9)A^a p U{9)+iU ] {9)d l U{9)) 

for £j. Thus £*,A? and IT, 0n are the ten pairs of physical variables. It therefore 
makes sense to choose £ and II as generators for the second-order Langevin system. In 
this case, however, a convenient system of equations is obtained by using a nontrivial 
multiplier Ili?^ 7 for the generators With this choice terms proportional to the noise 
take especially simple form 

ut]r^{£?, A?} = nr° (17) 

for the gauge potential. For the scalar field the corresponding term is 

nrj^{^,0} + r n {n,0}, (is) 

from which £, R, IT can be readily eliminated in favor of the original variables by noting 
that an arbitrary variation of the fields can be written as 

SAf = nr°, 50 = (p + iZ a a a )7i* (19) 



with real T,p, E, provided \n\ > (generic case). Gauss' law (|T5| ) is only respected if 

S a = -e a/37 lf El (20) 

Since variations ( ff7^l8 ) also respect the constraints, they must be of the form ( |19| , pOD . 
Comparing variation of 0n due to (|19D with that due to ( pT8[) yields p = T n . Finally, 
equality of terms linear in Tf gives 

UR^{£l<P} = -ie a ^E2(j a Tx*. 
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4 Lattice Langevin equations: formulation 



The next step in our program is space discretization of constraint-preserving thermal- 
ization scheme. Langevin equations of motion for lattice gauge theories are constructed, 
following the analogy with the continuum case. Having in mind applications to elec- 
troweak theory and to quark-gluon plasma, I will only write down lattice equations for 
the two nonabelian theories already considered. The reader should have no difficulty in 
writing a similar set of equations for the abelian theories. 

A convenient formalism for Hamiltonian (A = 0) lattice gauge theory is that of Kogut 
and Susskind |17]]. The configuration space consists of unitary matrices U\ forming the 



fundamental representation of the gauge group on every link / of a cubic lattice. In the 
following I shall use the notation j, n for a link along a positive direction n originating at 
site j. Lattice analogs of electric fields are link variables Ef** generating right covariant 
derivatives on the group. Poisson brackets obeyed by Ep 01 with U\ and among themselves 
are only nonzero for variables residing on the same link, and then, in the case of SU(2) 
gauge group 

{Ef a ,U l } = -iU l a a ; {E^, E^} = 2e a ^ E*\ (21) 

where a a are Pauli matrices. The link variables Ui and E^ a span the phase space. 
Alternatively, one could replace Ef** as independent variables by Ef 10 , the generators of 
left covariant derivatives on the group. On every link I the transformation between the 
two sets of variables reads 

E^ a a a = -E* P Uia p U}. (22) 

As a matter of convention, I choose here E^ a as independent variables. A useful property 
is 

{E* a , E^} = 0. (23) 

Scalar doublet fields <f>j, if added to the theory, are assigned, together with their conjugate 
momenta iTj, to the sites j of the lattice. Gauge transformations of the variables sharing 
site j are generated by 

°t = - £ [ e ja + E f-nA - 1 fao'h ~ ( 24 ) 

h 

(the scalar-field terms should be dropped from (^4|) in case of pure Yang-Mills theory). 
The set of Gauss' laws is CJ = 0. 

Consider a Langevin system for the Yang-Mills theory first. Proceeding by analogy 
with the continuum theory, I choose the gauge-invariant generators 

with 1 < n < n' < 3, and a unit multiplier. The value of the friction coefficient 7# > is 
arbitrary and can be tuned to optimize the algorithm performance. With the generators 
and the multiplier chosen, the Langevin system is fully determined by the choice of 
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a Hamiltonian. The SU(2) Yang-Mills theory on the lattice is usually described by a 
Kogut-Susskind Hamiltonian 

Hym = I E E^E^ + E (i " l^Uu) • (25) 

The second term in ( ^5|) is the standard plaquette term describing the magnetic part of 
the energy. The Langevin system can now be written explicitly using (pT0| . pT| , |22|^3| , |25|) 
(the site index j common to all variables is dropped to simplify the notation) 

Un = -tEtUna a -t^Y,{ Glhn ' ] + Tlhn ' ] ) E^^U h , (26) 

n' 

Here T n ,n' is an independent random variable corresponding to a given combination of 
a site j and [fin 1 ] , 

ee -!M^Tr (jg#V E + n <- n') ; G nn = *^^Tr (e%o« E , 

(27) 

and n h is any plaquette containing the link XJ<~ n . The E R equation of fl2TS|) has no stochastic 
or dissipative terms because these variables commute with the generators. Note that (f2~6[), 
beside satisfying the Gauss' law, automatically preserves the unitarity of link matrices, 
as it should. 

With the scalar field included, a Langevin system is again constructed by analogy 
with the continuum case. Gauge invariant generators n = |7r| 2 are introduced through 

The equations take a relatively simple form with multipliers ^/jUR^ for and y / 7n 
for n, where 7 and 7n are arbitrary positive friction coefficients. With an addidion of 
the scalar field a representative Hamiltonian is 

H H = H YM + E W 2 + E \^+n ~ UlM 2 + A E W (|0/) , (28) 

3 j,n j 

where W is a local scalar field potential. In the resulting Langevin system the ir and E R 
equations are simply the Hamiltonian ones, whereas for U and 4> one obtains 

Uj t h = 



Etn + v^Kf (r^,(t) + G&JJ t^ff"; 
* + v^rT (if + G°) tt* + E + G%) U^o*^, (29) 
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where, as before, T are mutually uncorrelated white noise variables, and 



G?fi = ^V7(^^ a + ^^ ft ^ 7 ^) +cc; 
Gf = -p^-^-^.+cc (30) 

(summation over SU(2) spinor and adjoint indices only). 

In the continuum theory, discussed in Section 3, the generators commuted among 
themselves and therefore could be made, by a suitable canonical transformation, canon- 
ical physical momenta. As a result, the continuum Langevin systems were second-order 
systems in the physical subspace. Since the lattice generators do not commute, the lat- 
tice Langevin system is not necessarily second-order in this sense. Nevertheless, in the 
low-temperature regime (/? > 1) the continuum Langevin systems are recovered. Indeed, 
the lattice generators do not commute due to the second relation of (pip . However, at 
low temperature E R ~ P -1 ^ 2 , as can be seen from fl2o]), hence the Poisson bracket of 
two E R s is (9(/5 _1//2 ), much smaller than {p, q] = 1 for any pair of canonically conjugate 
p, q. Therefore, at low temperatures both ( p6|) and (p9|) reduce to second-order Langevin 
systems in the physical phase space. 



5 Lattice Langevin equations: numerical integration 

The final step in our construction of thermalization algorithms is finding a suitable 
numerical integration scheme for a system of stochastic differential equations of the type 
([10|). It is best to use a scheme which respects the Gauss' law. Such a scheme is 
indeed possible due to the following useful property shared by all the Langevin systems 
constructed here. Namely, with the conventional choice of variables (space components 
of the gauge potential and the electric fields) ithe Hamiltonian is a sum of a kinetic term 
K commuting with all the momentum variables p (of which the lattice color electric 
fields are a generalization), and a potential term V commuting with all the coordinates 
q, with K and V being each gauge-invariant. Moreover, all the proposed generators 
T k of ( 10 ) commute with all the momenta. As a result, the Gauss' constraints C a 



are conserved separately by the momentum and coordinate equations of motion. This 
property was used, although not explicitly stated, for purely canonical equations in 
Ref. @. Moreover, since for the noise variables T in (^,^) all realizations in time are 
allowed, it is clear that this property holds if G terms in P3J20) are arbitrary functions 



of time, not necessarily given by fl27 ) and by (0). We conclude that the Gauss' law 
holds exactly for any combination of (a) integrating exactly the E an 7r equation of 
motion while keeping U and fixed, and (b) integrating exactly the U and <fi equations 
while keeping E R , n, and G fixed (the random variables T are also kept fixed for the 
duration of step (b)). The freedom in combining (a) and (b) can be used to optimize 
accuracy and stability of the algorithm. 
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Some technical adjustments may have to be made in applying this general idea to a 



specific gauge theory. Consider first the equation (pq) for the SU(2) Yang-Mills theory. 
For arbitrary T and G, the stochastic terms on the r.h.s. involve in a nonlinear fashion 
all the link variables for links emanating in positive directions from site j. Performing 
step (b) for such a system of coupled nonlinear equations is a very complex task. If at 
any step we only allow a single off-diagonal component of the symmetric tensor r + G to 
be nonzero, a considerable simplification will follow. Note that the general form of 
is (no implicit summation over n and n' in the following) 

tin = iY Z~ ~,E L ~,U n , 

•<■ / j nn' n ' 71 > 



where Z nn ,, symmetric under exchange of n with n', is constant for the duration of step 
(b), and E L n = E La na a obeys for fixed E R an equation of motion 

E L n = iY Z~~, \E L h, E L ~, 

/ j nn' ' n 



so that Yjn E L h = 0. Consider now a special case where the Z matrix has only one 



Zfii Zy2 



Z, Z 
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J 2H 



0. Then EL and 



3 - 



nonzero off-diagonal element, e.g., Z nn 

equvalently, E^ + E~ are time- independent, and the equations for E L and for U are easy 
to solve: 



Unit) 



exp(-iZ(Ef + E$)t)U n (0) exp(i(Z - Z 1 - 2 )E^a a t) 



U^(0)exp(iZ^ a a a t). 



In performing step (b) for the system (|29|) it helps to note that the third term in 
the (p equation is completely fixed by requiring that it compensates the violation of the 
Gauss' law generated by the U equation. This suggests the following procedure. First 
the equations for the link matrices are integrated. Then violation <5C° of the Gauss' law 
at site j, resulting from this change of link matrices alone, is determined. Finally, the 
Gauss' law is restored by the change in the scalar field: 



<f>(t)-<f>M 



1 + /H(r? + Gf) + -^- 2 e a(3l 5Cy a 



(31) 



There is one caveat in this procedure: if \irj\ is close to zero, divisions by required for 
the scalar field updates in ( |3TD can be numerically unsafe. At the same time, comparison 

1 7Tj | 2 approaches a finite value as 



of ( pTD with the equation of (p9|) shows that SC. 



7Y.; 



7To 



0. The remedy is therefore to set a minimal value e of 



7T; 



and, whenever 



< e, replace SC"/\iTj\ by an estimate obtained in the following way. First 



explicitly appearing in the Uj jn equations of ( p9| ) is replaced by e. Next, the \J^ n equations 
are integrated in order to find SC®. The estimate is then SC^/e. It is not difficult to 
show that, as a result of this replacement, the Gauss' law is maintained approximately, its 
violation being at most 0(e 2 ), which can be made completely negligible by an appropriate 
choice of e. 
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Steps (a) and (b) must be combined in a suitable fashion to form an accurate and 
stable integration algorithm. A simple scheme accurate to the order 3/2 in the time step 
A is as follows. Initially (t = 0) a set of random variables T is generated. These remain 
constant for the duration of the step. The integration step is schematically represented 
as 

1. E R (0), tt(0) -> E R (f ) , 7T (f ) for fixed U(0), 0(0); 

2. £7(0), 0(0) - U (f ) , (f ) for fixed (f ) , vr (f ) , G(0); 

3. C/(0), 0(0) - [/(A), 0(A) for fixed E R (f ) , vr (f ) , G (f ) ; 

4. £ fi (|) , 7T (f ) -> £ K (A), tt(A) for fixed £/(A), 0(A). 

Here the — > symbol denotes exact integration. All the color and space indices have 
been dropped in order to simplify the notation. For a pure Yang-Mills theory the 
and 7r variables should be omitted. Note that this scheme is a generalization of the 
simplest leapfrog algorithm to which it reduces in the absence of stochastic terms (item 2 
above then becomes obsolete). If not for the V terms present in the stochastic equations 
the algorithm would have an 0(A 3 ) error per step, similar to the simplest leapfrog. 
However, the noise terms T must be implemented as random variables having average 
and variance 2/A in order to approximate (0). The error per step is therefore 0(A 3 ^ 2 ). 

Numerical tests, conducted at high (/3 = 2) and low (/? = 12) temperatures, confirm 
that the algorithms in question indeed achieve the correct thermalization while main- 
taining the Gauss' law with high accuracy. The tests were performed on 12 3 lattices 
with algorithm parameters je — 0.05, A = 0.01 for the Yang-Mills theory and 7 = 0.04, 
7n = 0.2, A = 0.005 in presence of the scalar field. The scalar potential was A(|0| 2 — v 2 ) 2 , 
i.e., Hh corresponded to the bosonic SU(2) sector of the Standard Model. The algorithm 
was tested with A = 0.5, setting unit Higgs to vector boson mass ratio, and v 2 = 0.05, 
corresponding to the Higgs inverse mass of about 4.5 lattice spacings. 

Figure [l] illustartes the thermalization process for the SU(2) Yang-Mills theory de- 
scribed by fl25|). The initial configuration had zero electric fields and the energy far below 
the average one for the assigned inverse temperature (3 = 12. The system then rapidly 
reached thermal equilibrium. An extremely small static color charge per site generated 
by the algorithm is entirely due to a limited computer accuracy and does not depend on 
the time step. On a Cray-C90 processor with single-precision arithmetic the magnitude 
of the spurious static charge was less than 4 x 10 -12 per site at the end of evolution 
shown in Figure [I|. Similar negligible amounts of constraint violation were observed in 
tests of the Langevin system ( p9|) with the scalar field. 

Several criteria were used to judge how close the algorithms come to generating 
canonical distributions for the two theories in question. In one instance, the thermal 
average of an observable is known analytically. Namely, the radial component of the 
scalar field momentum, 2Re7r0/|0| appears in H H quadratically and can be chosen as 
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Table 1: Summary of measured canonical averages for spatial Wilson loops of various 
sizes at the beginning (initial) and at the end (final) of a Hamiltonian trajectory 50 time 
units long. 



an independent physical canonical variable. Hence the average radial kinetic energy per 
lattice site should be 1/2/3. This exact result was indeed found within margins of the 
measurement error in both low and high temperature regimes, with a 0.2% accuracy. 

Unfortunately, I am not aware of other observables with exactly known thermal 
averages. For some quantities, however, perturbative estimates should be reliable at low 
temperatures, opening another possibility of comparison with numerical experiment. In 
particular, since the low-temperature system is only weakly nonlinear, the average energy 
per degree of freedom should be close to the equipartition value 1/(3. At (3 — 12 the 
measured average energy per degree of freedom is 1.0091(7)//? for the Yang- Mills theory 
and 1.0053(7)//? with the scalar fields. 

Finally, the equilibrium statistical weight, exp(—(3H), being a function of the energy 
only, is conserved by the Hamiltonian evolution. This property is used for the follow- 
ing consistency test of the algorithms. An initial thermal configuration is subject to 
the Hamiltonian evolution, and observables are measured at the beginning and at the 
end of the Hamiltonian trajectory. For any observable, the two sets of measurements 
should yield the same average if the generated distribution is indeed canonical. This 
test was conducted for a number of observables of the Yang-Mills theory, namely, the 
color electric energy density, the topological charge density, and spatial Wilson loops 
of various sizes. With the scalar field included, the radial kinetic energy density, the 
Bricmont-Frohlich correlation function [IT9J at various distances, and the "string bit" 



(j)jUj t n(j)j + n/ \ |0j + n|) 1//2 were added to this list. All the measurements show that the 
generated distributions are close to the canonical ones. As an example, Table [I] summa- 
rizes the results for the spatial Wilson loops. 
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Figure 1: Langevin time history of the color electric energy per site for the SU(2) Yang- 
Mills theory on a 12 3 lattice. The algorithm time step was 0.01. Initial configuration 
had energy small compared to the average one for the inverse temperature j3 = 12. 

6 Application: maximal Lyapunov exponent of the 
SU(2) Yang-Mills theory 

As pointed out in the Introduction, lattice theories in the classical approximation are 
especially valuable in dealing with dynamical properties at finite temperature since in this 
case the arsenal of nonperturbative tools applicable to a full quantum theory is extremely 
scarce. Lyapunov exponents of a gauge theory belong to this category. According to 
recent numerical work ||, the maximal Lyapunov exponents A max of the classical SU(2) 
and SU(3) lattice gauge theories are approximately independent of the lattice spacing 
and proportional to the average energy per degree of freedom in a range of values of the 
latter. This approximate scaling property of A max was suggested based on a set of initial 
configurations with zero electric energy and the values of link matrices randomly chosen 
from a closed vicinity of identity. Here I use the thermalization algorithm for the SU(2) 
Yang-Mills theory described by fl25|) to compute (A max ), the canonical ensemble average 
of A max (notation () is used for canonical averages in the following). Taking this average 
allows to directly establish a relation between the maximal Lyapunov exponent and the 
temperature and to study deviations of (A max ) from the scaling prediction of Ref. || as 
compared to the statistical error of (A max ). 

The procedure used for measuring A max was as follows. A sample of phase-space 
configurations corresponding to an inverse temperature f3 was generated. Each member 
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10 20 30 40 50 60 70 80 

time 



Figure 2: Distance Dm between two diverging trajectories as a function of time for 
(3 = 2.5 (diamonds), (3 = 4 (squares), and (3 = 7.5 (pluses) on a 12 3 lattice. The error 
bars are smaller than the plotting symbols 

of the sample served as a reference initial configuration. A neighboring configuration 
was then generated by applying the algorithm briefly (in terms of Langevin time) to a 
reference one. Every such pair of nearby initial configurations was let evolve according to 
the Hamiltonian equations of motion, and the distances between the two configurations, 
defined as 

D E = J2 IW - E'?E'?\; D M = £ \TtU D - Tr^|, 
i □ 

||, were monitored (the unprimed and primed variables correspond to the reference 
and the neighboring configuration, respectively). A Langevin trajectory between two 
consecutive Hamiltonian ones was long enough to eliminate autocorrelations of (De^m^)) 
within the sample. The measurements were performed for inverse temperatures 2 < (3 < 
10. A 20 3 lattice was used for /3 > 5, and a 12 3 lattice for higher temperatures; this 
size reduction had no measurable effect on (A max ) already at (3 = 4. The Hamiltonian 
time step was small enough to ensure conservation of energy to six significant digits. For 
every value of the temperature 25 Hamiltonian trajectories were performed. 

As functions of time, {D e ,m{^)) both exhibit transient effects early on, then grow 
exponentially, and, finally, saturate (Figure |2|). The maximal Lyapunov exponent A max 
is defined as the rate of the exponential growth of De,m{P)- A fit of I?s,m(0 to an 
exponential function of time yields (A max ). Within the measurement errors the values of 
(A ma x), determined from both definitions of distance, coincide. Measurements in Ref. 
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{H YM )/V 



3.5 



Figure 3: Maximal Lyapunov exponent A max as a function of average energy per site. 
Diamonds with the error bars smaller than plotting symbols are canonical averages of 
the present work. Pluses and squares are single-trajectory results of Ref. [6]. The dotted 
line is a linear fit through the origin. 
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Table 2: Summary of measured canonical averages of A 



suggested scaling of (A max ) with the average energy per site: 

(Amax) / {Hym) = k. (32) 

According to Ref. k ~ 2/9 in the units of energy and time used in this work. The best 
linear fit through the origin for the data in Figure || is indeed very close to 2/9. However, 
the goodness of fit is very poor since, as Table shows, statistically significant deviations 



from fl32|) are as large as 10%. Hence the estimate of k in the range of temperatures 
considered may not be reliable. It is obvious from dimensional considerations that any 
deviation from scaling (A max ) oc T implies dependence of (A max ) on the lattice spacing. 
Lattice artifacts tend to increase with the temperature due to the compact form of the 



lattice magnetic term (cf (|25|)). Thus, in order to safely establish the value of k, A max 
should be measured at temperatures T < 0.1. Larger lattice sizes may be required 
in order to control finite size effects, clearly visible already at T = 0.1. An accurate 
knowledge of k through a low-temperature measurement is necessary in order to test the 
relation of (A max ) to the plasmon damping rate as suggested in Ref. jTj. 



7 Concluding remarks 

Methods presented in this work enable one to accurately bring a variety of classical lat- 
tice gauge theories to thermal equilibrium while respecting the local charge conservation. 
Due to their local nature, these methods easily lend themselves to vector and parallel 
computer implementations^. The algorithms can be applied to the study of real-time 
properties of high-temperature gauge theories, where the classical approach is by far 
the simplest if not the only one available. There is a growing body of evidence that 
some interesting nonperturbative quantities can be reliably determined by classical real- 
time simulations and are not plagued by ultraviolet divergencies inherent in classical 
thermodynamics. The examples include approximate scaling of thermalization rate in 

1 The codes are available upon request from the author 
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(3+l)-dimensional Yang-Mills theories discussed in the previous section, as well as nu- 
merical evidence for the existence of continuum limit of fermion-number violation rate in 
(l+l)-dimensional Abelian Higgs model |2], ||, |]]. More recently, the SU(2) Yang-Mills 
thermalization algorithm described here was applied to determine the baryon-number 
violation rate in the Standard model at temperatures well above the electroweak phase 
transition [|J. The results show the existence of continuum limit for the rate and its 
sensitivity to long-range properties of the theory 

It is possible to generalize the method to other gauge theories, most notably to the 
physically important SU(3) Yang-Mills theory. Again, the algorithm construction begins 
with finding a suitable set of gauge-invariant generators. A candidate generator set is 

= E a E a - = d 01 ^ 1 E Q E^ ''E 1 

i j i i j k ' 

where E are the color electric fields, d are the symmetric SU(3) constants, and i < j < k. 
Gauge invariance of T can be explicitly verified, they commute with each other, and their 
number is 16, the number of the physical degrees of freedom. A lattice algorithm can be 
developed along these lines. 

Finally, it should be noted that this work has concentrated on constructing the al- 
gorithms and verifying their validity, rather than on tuning the algorithm performance. 
The latter is postponed to subsequent work. For the most obvious and important applica- 
tions of the algorithms, however, namely, real-time studies of gauge theories, generation 
of initial thermal configurations typically consumes little time compared to canonical 
real-time evolution. With this kind of application in mind, the algorithm optimization 
can only lead to modest gains in computing time. 
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